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Abstract 

We study the traffic of two types of molecular motors using the two-species asymmetric simple 
exclusion process (ASEP) with periodic boundary conditions and with attachment and detach- 
ment of particles. We determine characteristic properties such as motor densities and currents by 
simulations and analytical calculations. For motors with different unbinding probabilities, mean 
field theory gives the correct bound density and total current of the motors, as shown by numerical 
simulations. For motors differing in their stepping probabilities, the particle-hole symmetry of 
the current-density relationship is broken and mean field theory fails drastically. The total motor 
current exhibits exponential finite-size scaling, which we use to extrapolate the total current to the 
thermodynamic limit. Finally, we also study the motion of a single motor in the background of 
many non-moving motors. 
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I. INTRODUCTION 



Intracellular transport by molecular motor proteins is essential for many cellular functions 
such as shuttling vesicles or organelles to their proper destinations [l, 2]. In eukaryotic cells, 
polar cytoskeletal filaments serve as tracks for three families of cytoskeletal motors: kinesins 
and dyneins that move along microtubules, and myosins that walk along actin filaments 

13. 

The movements of a molecular motor along its filamentous track can be described as 
a biased random walk. However, after a certain number of steps the motor unbinds from 
the track because the thermal noise of the environment overcomes the finite motor-filament 
binding energy. The motor then diffuses in the surrounding fluid as a Brownian particle 
until it comes close to the track to rebind to it. On large time scales, motors therefore 
perform peculiar 'motor walks' that consist of alternating periods of directed motion or 
biased random walks on the tracks and undirected diffusive motion in the fluid environment 

a. 

Intracellular transport is typically achieved by several types of motors rather than just 
by one type . Different types of motors may either move independently along the same 
filament or work together to power bi-directional transport of a cargo that carries two types 
of motors [6( . Here we consider the first case, where different species of motors move on the 
same filament, and interact only via their mutual exclusion. These motors species might 
move into different directions, or they might move into the same direction but have different 
velocities or different affinities to the filament. We are interested in the traffic behavior of 
the latter systems. 

A simple model for motor traffic is based on the asymmetric simple exclusion process 
(ASEP) [7]. In the ASEP, particles that interact through mutual exclusion from lattice sites 
perform biased random walks along a one dimensional lattice. This paradigmatic model for 
non-equilibrium transport exhibits interesting phenomena such as boundary induced phase 
transitions and traffic jams [7]. To describe the traffic of motor proteins, ASEP-like models 
have to be supplemented with the dynamics of binding to and unbinding from the track to 

account for the finite run length of the motors flat' 
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IPs with a single species of particles can be generalized to models with multiple species 



13lll4|. Even without attachment and detachment of particles, multi-species ASEPs exhibit 
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a variety of cooperative phenomena, such as spontaneous symmetry breaking 15J and phase 



separation 
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To model cellular traffic with different types of motors, one has to consider multi-species 
ASEPs with attachment and detachment of particles. Apart from investigations on condi- 
tions for the exact solvability jisl . , most work on multi-species ASEPs with detachment 



and attachment has 
directions 
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ocused on the traffic of two species of particles moving into 
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rese models exhibit continuous phase transitions 



hysteresis of the total current 
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21] . traffic jams 24|, spontaneous symmetry breaking and 



localized shocks 221 ] . Here we focus on the case where two species of particles move in the 
same direction, but do so either with different velocities or with different binding affinities 
to the filament. A few such systems have been studied experimentally in vitro, in particular 
the traffic of mixtures of fast-moving and slow- or non-moving kinesin motors [25]. 

The paper is organized as follows: In section [Til we will present the theoretical model for 
multi-species motor traffic. We will focus on the situation that two motor species differ in 
only one parameter. In section II III we will discuss the traffic of two types of motors which 
differ only in their unbinding parameters. In section IIVI we will discuss the traffic of two 
species with different stepping parameters, specifically one species with zero velocity. We 
will end with a summary and discussion in section [V] 



II. THEORETICAL MODEL FOR MULTI-SPECIES MOTOR TRAFFIC 

To investigate cellular transport or in vitro experiments with multiple species of molecular 
motors, we study multi-species ASEP models with attachment and detachment of particles. 
As mentioned in the introduction section, motors perform random walks, which are biased 
towards the same direction for all motor species, along a cytoskeletal filament. Such filaments 
are polymers which have periodic motor binding sites with repeat distance £, and motors 
walk along the filament with steps of the same size I. For microtubules, the lattice constant 
is t = 8 nm. We therefore represent the filament by a one-dimensional lattice of L binding 
sites with spacing i. Per unit time r, a bound motor of type % on the filament makes a 
forward step with probability a\ if the target site is free, unbinds with probability e; and 
remains at the same site otherwise [8(. Unbound motors are treated as a motor reservoir 
jsl, [h| with unbound density p u b,i for the z-th type of motors. An unbound motor binds to 
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the filament with the sticking probability 7r a d,i when the binding site is empty. The treatment 
of unbound motors as a reservoir is slightly different from that in a previous lattice model for 



motor traffic in Ref. [8[] where the diffusion of unbound motors has been treated explicitly 
as Brownian movement. The simplification to treat the unbound motors as a reservoir is 
appropriate if the density of unbound motors is homogeneous, as is the case for sufficiently 
fast diffusion and large motor numbers in the bulk solution. This is typically the case in in 
vitro experiments, where the buffer is an aqueous solution with motor concentrations in the 
nano- or micro-molar range. The typical length of a microtubule is of the order of 10/im 
or even larger, which is one order of magnitude larger than the typical motor run length, 
which is of the order of 1/xm for kinesins. Therefore, we use periodic boundary conditions 
throughout this paper. This simplification is appropriate as most motors do not feel the ends 
of the filament, so that effects of the filament ends can be neglected. 1 In our simulations, 
the length L of the filaments is chosen sufficiently large so that finite-size effects can be 
ignored, as we show explicitly in Sec. [IV] below. 2 

% PjbJ 

— • — • — • — • — • — • — • — 

FIG. 1: (Color) The model for multi-species motor traffic. Per unit time r, bound motors of species 
i hop one step forward along the filament with probability a; and detach from the filament with 
probability ei, while unbound motors with density p u h,i attach to an empty site with probability 
7r a d ; i. Exclusion interactions between motors prohibit the movement of motors if the target site is 
occupied by another motor. The system is on a ring of L sites with lattice constant t. 

To sum up, as shown in Fig. [TJ, during each time step r one of the L sites is randomly 



1 Note also that we do not scale the binding and unbinding probabilities of motors with the system size L. 
The latter scaling, which has been used in Ref. [n]], leads to a motor run length that is comparable to 
the filament length L in the limit of large L. 

2 In most of our simulations, the filament length L was chosen to be L = 200 corresponding to about 2/im. 
In the presence of slow motors, the actual run length of the fast motors is strongly reduced and small 
compared to L = 200. In the absence of slow motors, as in Sec. IIII1 the length of the filament is irrelevant 
as well since our system then exhibits mean-field behavior for all L. 
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chosen and updated according to the following stochastic dynamical rules: 

Stepping: iO— >0i with probability (1) 
Detachment: i — > with probability e\, (2) 
Attachment: — > i with probability p u b,i7i"ad,i, (3) 

where denotes a vacancy site and % represents a site occupied by a motor particle of species 
i. All probabilities for other transitions are zero. 

We performed Monte Carlo simulations on the the one-dimensional periodic lattice with 
L sites. In the simulations, we control the number of motors in the system by changing the 
value of the unbound motor densities. This mimics in vitro experiments where the motor 
concentration in the buffer is an easily accessible control parameter. Cells can also control 
the motor concentration in the cytoplasm by regulating gene expression. Once the unbound 
motor densities have been settled, we use a random-sequential update. Each simulation step, 
which corresponds to a unit of the basic time scale r, consists of L Monte Carlo moves. At 
each such move, a lattice site is chosen randomly and updated according to the dynamical 
rules ([T]) to The simulations are performed long enough (10 10 steps) to ensure that the 
steady state of the system has been reached. The current and bound density are calculated 
by averaging over all the simulation steps and over all lattice sites. 

In this article, we consider the traffic of two species of motors for long times, i.e. in the 
stationary state. Generally speaking, all the parameters for different motor species might 
be different. Here, as a simple starting point of theoretical work, we focus on the situation 
that the two motor species differ in only one parameter. We consider different unbinding 
probabilities, i.e. e\ ^ €2, in section [TTT1 and different forward stepping probabilities, i.e. a± 7^ 
«2, in section HVl We will not explicitly consider the case of different binding probabilities 
for different motor species, i.e. 7r ad l 7^ vr ad 2 , because this case is equivalent to that of 
different species with only different unbound densities, since the attachment probability is 
the product of binding probability for single motor 7r a d,i and unbound motor density p u b,i, 
see the dynamical rule fl3]). 

Since we are interested in modeling molecular motor traffic, we choose the probabilities 
for motor movement as found experimentally for some specific molecular motors. We focus 
on the well-studied motor kinesin-1 for which all necessary parameters have been measured 
in single molecule experiments. Kinesin takes about 100 steps before unbinding from the 
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TABLE I: Summary of the typical values of probabilities per time step r = 10 4 s used in this 
article, which mimic the properties of the molecular motor kinesin-1. 
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FIG. 2: (Color online) Traffic of two motor species with different unbinding probabilities, (a) 
Total bound density and (b) normalized total current J /a (in the unit of t -1 ) as a function 
of total unbound density p u b for different fractions of type-2 motors. The simulation results 

(data points) agree well with the mean field calculation of Eqs. (JHJ), ([9]) and (|10p (solid lines). The 
parameters are e\ = 10~ 4 , €2 = 8e%, a = 10 -2 , 7r a d = 1 and L = 200. 

filament jlj , so that we can estimate the ratio of stepping and unbinding probabilities to be 
ale ~ 100. Kinesin's microtubule desorption constant -Kd,mt is of the order of 0.1 — 1 fiM 25 . 



Dinding probabilities e/7r ad = Na^ 3 Kd >mt 



271 ] . We choose the discretization time unit 



261 ]. which fixes the ratio of the unbinding and 
10~ 4 with the Avogadro constant A^a, see Ref. 
t so that the binding probability 7r a d = 1, which means that we have unbinding probability 
e = 10~ 4 and stepping probability a = 10~ 2 . Considering kinesin's velocity = aljr w 
1/xm/s [l|, this means that our discretization time is r pa 10~ 4 s, which is sufficiently small 
to avoid discretization artifacts. The typical values of the probabilities used in this article 
are summarized in Table [B 
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III. TWO SPECIES OF MOTORS WITH DIFFERENT UNBINDING PROBA- 
BILITIES 



In this section we consider the case of two species of motors which have different unbinding 
probabilities e\ and eg. The stepping probability a and the binding probability 7r a< j for both 
species are the same so that the indices of the corresponding symbols can be omitted. 

Since the detachment and attachment rules ([2]) and ([3]) do not involve neighboring sites, 
there is no correlation of motors in the processes of binding and unbinding. Therefore the 
balances of attachment and detachment for both species of motors in the stationary state 
can be described by mean field theory as follows: 

VadPub,^ 1 - Pb) = ClPb.l, (4) 
7T a dPub,2(l - Pb) = e 2 Pb,2, (5) 

where pb,i and pb,2 are the bound densities of type-1 and type-2 motors, respectively, and 
Pb = Pb,i + Pb,2 is the total bound density. The terms on the left-hand side of Eqs. 01]) and 
(JI]) represent the attachment of motors, with the mean field factor (1 — pb) for the exclusion 
interaction between motors. The terms on the right-hand side represent the detachment of 
bound motors. Here we ignore the exclusion interaction when motors try to detach from the 
filament since the total unbound density of motors p u b = p u b,i + Pub,2 is very low, typically 
below 10~ 3 , which corresponds to a bulk concentration of unbound motor proteins below 
Pu b /£ 3 N A « 1 M. 

For convenience, we define the fractions 

* ub ,i = ^ and * b)i = ^ (6) 

Pub Pb 

of the z-th type of motors in the unbound and bound states, respectively, with ^ub.i + ^ub^ — 
^b,i + ^b,2 = 1- Because the two species of motors have different unbinding probabilities, 
the fractions of these two types of motors in the bound states are not equal to those in the 
unbound states, or more precisely, 

*b,l £2 ^b,2 ^ _^2_ ^ 



*ub,l Cl ^ub,2 #ub,2 

From equations fll]) and ([SD we obtain the bound densities of both types of motors and 
the total bound density 

^adPub fQ s 

Pb = (8) 

TTadPub + ^eff 
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with an effective unbinding probability 

£1 £9 

e eff = ei^b,i + e 2 ^b,2 = — Tf ■ ^ • (9) 

Eq. (|S]) has the same form as the relation of binding and unbinding probabilities in the 
1-species ASEP with attachment and detachment, except that the unbinding probability of 
the single motor species is replaced by the weighted average effective unbinding probability 
e e fj. In the limit of \l/ u b,i = or \l/ u b,2 = 0, the effective unbinding probability becomes 
e e fr = £2 or e e fj = ex, respectively, and Eq. (jSJ) reduces to the corresponding relation of the 
single-species ASEP. 

The total current J of both species can be determined by mean field theory as the product 
of forward stepping probability a, total bound density pb of the two types of motors and 
density of vacancies 1 — Pb, 

rv 

J=-p b (l-Pb). (10) 

T 

The total motor current has a maximum at pb = 0.5. It decreases for bound densities larger 
than p b = 0.5, because overcrowding of bound motors leads to traffic jams. 

Simulation results of the ASEP with two species of particles with different unbinding 
probabilities are in excellent agreement with mean field theory, see the example with t 2 = 8ei 
in Fig. [2J This observation strongly suggests that mean-field theory is exact in the case of 
two species of motors which differ only in the unbinding probabilities. 

Results of the mean-field theory for the ASEP with two species of particles with different 
unbinding probabilities are in excellent agreement with simulation results (see the example 
with €2 = 8ei in Fig. and with the numerical solution of the associated stochastic process 
for small system sizes (not shown). This observation indicates that mean-field theory is 
exact in the case of two species of motors which differ only in the unbinding probabilities. 

As shown in Fig. [21 the total bound density pb increases with the total unbound density 
p u b, as to be expected. The total current J as function of unbound density p u b has a 
maximum because overcrowding by motors causes traffic jams. The maximum shifts to 
higher unbound density with increasing fraction of type-2 motors, because type-2 motors 
have a larger unbinding probability than type-1 motors. Therefore a higher fraction of type-2 
motors leads to smaller bound density and less traffic jam. 
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IV. TWO SPECIES OF MOTORS WITH DIFFERENT STEPPING PROBABILI- 
TIES 



We next consider two motor species that differ only in the forward stepping probabilities, 
i.e. ai 7^ a 2 - The binding probability 7r a< j and unbinding probability e for these two species 
are the same so that the indices of the corresponding symbols can be omitted. Without loss 
of generality, we assume that type-1 motors have a larger stepping probability than type-2 
motors. We will focus on the case that the second motor species has zero forward stepping 
probability ct 2 = 0, because we expect that this largest difference between a,\ and a 2 leads to 
the strongest effect. This choice is also motivated by in vitro experiments, which addressed 



the traffic of mixtures of moving motors and non-moving mutant motors 251 ] . 



A. Traffic of two species of motors with different stepping probabilities 

As already discussed in section 1111} there are no correlations of motors in the processes 
of attachment and detachment. Therefore the balances of binding and unbinding processes 
for both kinds of motors can be described by mean field theory as given in Eqs. (jl]) and (J5]) 
with the simplified condition e x = e 2 = e. From these equations, the total bound density p b 
of the two types of motors is given by 

Pb = — ■ (11) 

TTadPub + C 

Eq. (1TT1) can be viewed as a special case of Eq. (JSj) with e\ = e 2 = e, so that e e fr = e, see Eq. 

Since type-1 and type-2 motors have the same binding and unbinding probabilities 7r ac j 
and e, the fractions of both types of motors in the bound and unbound states are equal, 

^ub,i = #b,i = % fori = 1,2, (12) 

with + * 2 = 1- 

The total current is given by the sum of the currents of the two motor species as J = 
Ji + J 2 . Within the mean field approximation, the currents of the two motor species are 

4f,i = — Pb,i(l - Pb) for i = 1, 2. (13) 
r 
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Note that in this approximation, both species experience the exclusion interaction through 
the same hindrance factor 1 — pb, which depends on the total bound motor density p^ So 
the mean field total current is 

7^7 7 - + ^2«2 n v 

%F — JmF,1 + ^MF,2 — PbU ~ Pb)- U 4 J 

r 

As a function of the total bound density pb, ^mf has a maximum at pb = 0.5 and is 
symmetric with respect to pb = 0.5. However, simulations show that the actual value of the 
total motor current for large systems is much lower than the mean-field prediction, see Fig. 
|3j In addition, the total current as function of bound density is no longer symmetric with 
respect to pb = 0.5 for L > 2. These deviations from the mean-field prediction indicate 
strong correlations of bound motors which come from the blocking of type-1 motors by the 
non-moving type-2 motors. 

To take these correlations into account, we consider the exact time evolution of the 
system. The possible configurations of the system can be labeled by the sequence of lattice 
site occupancies. For example, the configuration C = 012 of a system of size L = 3 means 
that the first site is empty and the second and the third sites are occupied by a type-1 motor 
and a type-2 motor, respectively. At time t, the system is in configuration C with probability 
P(C,t). At the next time step t + r, the system evolves stochastically following the rules 
(CQ), (T5]) and ([3]), and the probability P(C,t) to find the system in configuration C satisfies 
the Master equation 

P(C, t + r)- P(C, t) = Yl M ( C > C WC t)- (15) 

c 

Here, the off-diagonal terms M(C, C) of the transition matrix M with C ^ C represent the 
probability of a transition from configuration C to configuration C during the time interval 
r, and the diagonal terms 

M{C ) C) = ~Y. M{C '^) ( 16 ) 

represent the total exit probability from the corresponding configuration C. 

For long times, the system evolves into the stationary state with time-independent prob- 
abilities P st (C), which satisfy 

= ^M(C,C')P st (C), (17) 
c 
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FIG. 3: (Color online) Traffic of two species of motors with different stepping probabilities. The 
normalized total motor current J/ot\ (in the unit of t _1 ) as a function of total bound density p^. 
Mean field theory as given by Eq. (|14h (solid black line) predicts a much larger total current than 
observed by simulation of a large system of size L = 200 (data points). The other lines show the 
exact analytical results from solving Eq. (|17p for small system sizes L = 2, 3, 4, 7, with the result 
for L = 2 being identical to the mean field result. The parameters are ol\ = 10 -2 , «2 = 0, 7r a( j = 1, 
e = 10~ 4 and <£ 2 = 0.3. 

and are therefore given by the null eigenvector of the the transition matrix M. From the 
stationary distribution P st , it is straightforward to obtain the total bound density and total 
current by averaging over the appropriate configurations. 

For a system of size L with 2 motor species, there are in total 3 L possible configurations. 
Therefore the dimension of the transition matrix increases exponentially with the system 
size L and Eq. (fTTj) can in practice only be solved for very small systems. 

The dimension of the transition matrix can be reduced based on the symmetry of the 
system. Since the system is homogeneous with periodic boundary conditions, configurations 
which differ by a translational shift are equivalent. For example, the configurations 102, 
021 and 210 of a system of size L = 3 are all equivalent and can be combined into one 
equivalence class. Since the translation operators form a cyclic group Cl, the number of 
equivalence classes for arbitrary size L is equal to the cycle index Z(Ci), which is given by 
Burnside's Lemma 28] : 



where GCD(fc, L) is the greatest common divisor of k and L and the number 3 comes from 




(18) 



fc=i 
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FIG. 4: (Color) Traffic of two species of motors with different stepping probabilities. The normal- 
ized total motor current J/a\ (in the unit of r _1 ) decreases with the system size L to a plateau 
value Jqq. The simulation data (points) are well fitted by the exponential function of Eq. (|19p (solid 
lines) , which itself is a fit to the exact analytical solution of Eq. (|17|) for small systems (white x) . 
The inset shows a semi-logarithmic plot of the normalized total current J* ja\ = ( J — Joc)/ot± (in 
the unit of t -1 ) from which the plateau value has been subtracted. The parameters are a± = 10~ 2 , 
"2 = 0, 7r ad = 1, e = 10" 4 and f 2 = 0.3. 

the three possible configurations 0, 1 and 2 for every filament site. This symmetry allows 
a reduction of the matrix size by combing all equivalent configurations into corresponding 
equivalence classes and appropriately modifying the transition probabilities, as described in 
detail in the appendix. In the limit of large L, this procedure reduces the dimension 3 L of 
the transition matrix by a factor of 1/L, as can be seen from Eq. (j!8p . 

We then solve for the null eigenvector of the reduced transition matrix and obtain the 
exact total current for system of size L = 2toL = 10. As shown in Fig. [3l the total current 
for system size L = 2 is identical to the mean field total current (1141) as can be checked 
explicitly by exact solution of Eq. (TT7I) for L = 2. For systems with L > 2, the particle-hole 
symmetry on the filament is lost, and the asymmetry of the total current as function of the 
total bound density with respect to pb = 0.5 increases for increasing L, see Fig. [3j 

As can also be seen in Fig. [31 as well as in Fig. HI the total current decreases for increasing 
system size L. The reason is that in systems which are small compared to the motor run 
length, a stepping motor 'sees' the hole left by its image in front, which leads effectively to 
reduced congestion and larger total current. As the system size increases, this effect becomes 
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FIG. 5: (Color online) Traffic of two species of motors with different stepping probabilities. Nor- 
malized total motor current J /a\ (in the unit of r _1 ) as a function of total bound density pb for 
different fraction \&2 of type-2 motors. The extrapolations (lines) following Eq. (|19p from exact 
results for small systems of sizes L = 2 — 8 are in fair agreement with the simulation results (data 
points). The parameters are a\ = 10~ 2 , «2 = 0, 7r a d = 1, e = 10 -4 and L = 200. 

less pronounced, so that the total current decreases and reaches a plateau value for system 
sizes comparable to the run length, which is 100 sites for our choice of parameters, see Fig. 
H 

As shown in the inset of Fig. HJ a more detailed look at this finite-size effect reveals that 
the approach of the total current towards its plateau value can be approximated by the 
function 

J(L) = Joo + ba- L . (19) 

Here, the two parameters a and b, which characterize the exponential decay, depend on the 
dynamic parameters of the system such as forward stepping probability a,\ or the fraction 
\&2 of type-2 motors. 

We now use this finite-size effect to obtain the total current for large systems. We first 
obtain the exact results for the total current of small systems with size L=2 to 8 by solving 
Eq. ffTTl) . Then we fit these total current values by Eq. (fT9l to estimate the plateau total 
current J^. This extrapolation method gives fair agreement with the simulation results, see 
Figs. H and El 

Fig. [5] shows the results for the total motor current as function of total bound density. 
Since non-moving type-2 motors tend to hinder the stepping events of type-1 motors and 
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hence enhance the traffic jams, the current- density curves are asymmetric with respect to 
Pb = 0.5, with the position of the maximum shifted towards the low density region. As the 
fraction of type-2 motors increases, the total motor current decreases. This effect is quite 
dramatic: For instance, the total motor current drops more than one half as the fraction 
of type-2 motors increases from to 0.1. We note that a similar effect has recently been 
discussed in the context of the traffic of RNA polymerases in transcription of bacterial 
ribosomal RNA [2j|. 

Our finding of the exponential finite-size scaling of the total current is different from the 
1 / L-scaling observed in the single-species ASEP with periodic boundary conditions, which 



arises from the constraint of a fixed particle number in such a model 30]. The presence 
of attachment and detachment kinetics removes this constraint and leads to weaker finite- 
size effects. Accordingly, we do not observe any finite-size corrections for the case ^2 = 0, 
which corresponds to a single species ASEP with periodic boundary conditions and particle 
attachment and detachment (data not shown). Exponential finite-size corrections to the 
total current, as observed for \?2 > 0, are also obtained in single-species ASEPs with open 
boundaries in the low or high density phase [14|. It appears that the exponential scaling 
can be obtained when there are interfaces between regions of different densities or boundary 
layers, i.e. interfaces at the boundary of the system. This is corroborated by the fact that the 
finite-size effect disappears for the single species ASEP with open boundary conditions with 
or without particle attachment and detachment kinetics, if the boundary rates are chosen so 
that there are no boundary layers. It is likely that the local traffic jams induced by the slow 
particles for ^2 > have an effect similar to that of jamming at interfaces between regions 
of different densities. 

B. Single moving motor with obstacle motors 

From a theoretical point of view, single tagged particles have been used to study inter- 



esting phenomena in ASEP- type models, such as shocks [3JJ. The motion of a single labeled 
motor in the environment of many other motors has also been studied in single molecule 
experiments j25[ | . 

These experiments motivate us to consider the case of only one moving type-1 motor which 
has a nonzero stepping probability a\ in the background of type-2 motors with stepping 
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FIG. 6: (Color) Velocity i> (in the unit of £/r) of a single moving motor as a function of the 
unbound density p u b.2 of obstacle motors. The simulation results (data points) are in excellent 
agreement with the analytic results given by Eq. (|23p (lines). The parameters are ai = 0, 7r a( j = 1, 
e = 1(T 4 and L = 100. 

probability a<i = which we will from now on call obstacle motors. Both types of motors 
have the same unbinding probability e and binding probability 7r a( j. 

Because there is only a single type-1 moving motor, the total bound density given by Eq. 
f TTTj) is basically identical to the bound density of type-2 obstacle motors pt>,2 provided that 
the system is sufficiently large. 

Since we are interested in the movement of the bound type-1 motor, we will only take into 
account the bound state of the moving motor. By noticing that there is only one moving 
bound motor, we can label the states of the system by the site which is just in front of the 
moving bound motor. There are two states of this site: occupied by an obstacle motor or 
empty, and these states are denoted by 12 and 10, respectively. 

The system changes from state 12 to state 10 with probability e when the obstacle motor 
in front of the moving motor detaches from the filament: 

12 10 (20) 

Likewise, the state of the system changes from 10 to 12 with probability 7r a( iPub,2 when an 
obstacle motor binds to the empty site in front of the moving motor, and with probability 
Q!iPb,2 when the moving motor moves forward (with probability a\) and meets an obstacle 
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motor (with probability Pb,2}'- 

10 12 or 10X ^4 012 

Here X denotes a site that is occupied by an obstacle motor with probability pt>,2 or free 
with probability (1 — p b ,2)- 
The transition matrix 

-e 7r ad p ub ,2 + «lPb,2 
e -^adPub,2 - «lPb,2 



(21) 



for the state vector (12, 10) T determines the stationary state of the system as the null 
eigenvector 

E =( ^b.a + aiPb.a e \ T 

\ e + 7TadPub,2 + «lPb,2 C + 7T ad Pub,2 + "lPb,2 / 

Since only the state 10 contributes to the motion of the moving bound motor, its velocity 
is given as 

v=(o i a 1 i)E = ai ] m - - r . (23) 

As shown in Fig. [BJ the velocity of the moving bound motor decreases as the concentration 
of obstacle motors increases, and the exact Eq. (1231) is in good agreement with the simulation 
data. 



V. SUMMARY AND DISCUSSION 



Motivated by cellular traffic, which is powered by several different species of molecular 
motors, we have studied the traffic of two motor species along the same track. For this 
purpose we have used a two-species ASEP model with periodic boundary conditions with 
attachment and detachment kinetics, in which the two motor species move into the same 
direction but have different velocities or unbinding rates. We have determined characteristic 
properties of the movements, such as motor density and total motor current, especially the 
current-density relationship, by computer simulations and by analytical calculations. As in 
the cell or in in vitro experiments, our control parameter is provided by the densities of 
unbound motors. 

We have considered two motor species that differ in only one property. If the two motor 
species only differ in the unbinding probabilities, increasing the fraction of motors with 
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larger unbinding probability leads to a decrease of the total bound density of motors. As a 
consequence, upon increasing the fraction of motors with larger unbinding probability, the 
total current decreases in the low density region and increases in the high density region, 
see Fig. |2Jb). For the case of two motor species with different velocities we have focused 
on the extreme case of non-zero and zero stepping probability for the two motor types. As 
to be expected, the total motor current decreases with increasing fraction of non-moving 
motors. Since particle-hole symmetry is broken for two motor types with different stepping 
parameters, the current- density relationship becomes asymmetric with respect to the bound 
density pb = 0.5. As shown in Fig. [3l its maximum shifts towards lower bound densities 
because non-moving motors hinder the stepping of moving motors and hence enhance traffic 
jams. 

To obtain these results, we have employed computer simulations and analytical calcu- 
lations. We have first used mean field theory to analyze our model. The attachment and 
detachment processes only involve single filament sites and therefore do not introduce new 
correlations into the model, so that mean field theory works well for the case of two species 
which differ only in their unbinding parameters. However, mean field theory fails to estimate 
the total current of two motor species with different stepping probabilities, as illustrated in 
Fig. [3j In this case, the exclusion interaction leads to strong interference of the motors 
moving at different velocities and therefore to correlations and failure of mean field theory. 
In order to obtain analytical results in this latter case, we go beyond mean field and use 
finite-size scaling of the total current to extrapolate the total current value of small systems 
to the thermodynamic limit of large system size. This is possible because the total current 
decreases exponentially with the system size, as shown in Fig. HI 

This work presents a simple model for traffic by multiple motor species, which can be 
regarded as a starting point towards the understanding of intracellular transport with dif- 
ferent kinds of molecular motors. Our model could be extended in many ways, such as by 
considering more than two different motor species which differ in multiple properties. Also, 
the model for the motors can be refined to include internal states of the stepping process 
or interactions between the motors which go beyond spatial exclusion. In the simple form 
presented here, our model gives a qualitative picture of the effects of traffic by several motor 
species. It makes predictions for the current- density relationship (Fig. [2] and Fig. [5]) and for 
the velocity of a single motor in the presence of immobile obstacles (Fig. [6]). It also suggests 
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that the movements of different species of motors, and therefore of different kinds of cargo, 
can be cross-regulated, so that the movement of one type of motor is controlled indirectly 
by the presence or absence of other types of motors or by changes in the parameters of these 
motors. As an example, in the case of two species with different stepping probabilities, a 
small fraction of slow moving motors can strongly suppress the movement of fast moving 
motors and lead to a significant reduction of the current. This may be useful for the regula- 
tion of transport of different types of cargo along the same filament, both for in vivo traffic 
and for in vitro biomimetic applications. 
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APPENDIX 

In this appendix, we explain in detail how we calculate the exact total current of the 
system studied in section IV A. We illustrate the method for a system of size L = 2. We 
recall that the type-1 motors have positive stepping probability a± and type-2 motors have 
zero stepping probability a<i = 0. Both species have the same unbinding probability e and 
the same binding probability 7r a d. 

For a system with size L = 2, there are 3 2 = 9 configurations: 00, 01, 02, 10, 11, 12, 20, 21 
and 22, where 0, 1 and 2 denote an empty site, a site occupied by a type-1 motor and a site 
occupied by a type-2 motor, respectively. We label these 9 configurations by k — 0, . . . , 8 
and combine them to the state vector (00, 01, 02, 10, 11, 12, 20, 21, 22) . The corresponding 
transition matrix M, as defined in Eqs. ffT5j) and (fT6j) . is given by 
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(A.l) 



V 

The diagonal terms M^k are abbreviations for M(C^ k \C^ k '), k = 0, . . . ,8, which represent 
the exit probability from configuration during time interval r. In order to conserve the 
probability, the sum of M^ and all other elements in the column k is zero, compare Eq. 
(HSJ). For instance, M , = -27r ad p ub; i - 27r ad p ub;2 and M 1A = -e - a x - 7r ad p ub> i - 7r ad p ubi2 . 

The stationary state of the system P st is given by the null eigenvector of the transition 
matrix M as presented in Eq. ( fl7l) . By averaging all the contributions to the total current 
from the stationary probability distribution, we get the total current: 



J = Tr(D • P st )/rL, 



(A.2) 



with a diagonal matrix D, which represents the contributions of each configuration to the 
total current. For L = 2, this matrix is given by 



D = diag{0, oti, ct 2 , oti, 0, 0, a 2 , 0, 0}. 



(A.3) 



The dimension of the transition matrix for the all 3 L configurations increases exponen- 
tially as the system size L increases. With the help of translation symmetry, the dimension 
of transition matrix can be reduced by a factor ~ 1/L as explained in Eq. (TlBl . see also 
Table HO 

The translation operators form a cyclic group Cl- The quotient set of the cyclic group C 2 
contains i?(C 2 ) = 6 equivalence classes, i.e. [00], [01], [02], [11], [12] and [22], which contain 
1, 2, 2, 1, 2, 1 equivalent configurations, respectively. Each equivalence class represents 
all configurations which differ by a translational shift. For instance, the equivalence class 
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TABLE II: The translation symmetry reduces the dimension of the transition matrix 



system size 
L 


dimension of original 
transition matrix 


dimension of reduced 
transition matrix 


2 


9 


6 


3 


27 


11 


4 


81 


24 


5 


243 


51 


6 


729 


130 


7 


2187 


315 


8 


6561 


834 


9 


19683 


2195 


10 


59049 


5934 


11 


177147 


16107 


12 


531441 


44368 



[01] represents not only the configuration 01 but also the configuration 10 which can be 
translated into 01 by the cyclic group C 2 . 

The corresponding reduced transition matrix is obtained by adding up all the transition 
probabilities within the same equivalence class. For L = 2, the reduced transition matrix 
for the vector of equivalence classes ([00], [01], [02], [11], [12], [22] ) T is 

M' = (A.4) 
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with 



K,i = 


-27T ad p ubi i - 27T ad p ubi2 , 


M' = 

2,2 


— 2c — ziiuA Aih i — 27r ar i AiK 9 


^3,3 = 


-2e - 27T ad p ub> i - 27T ad p ub>2 


Ka = 


-2e, 


= 


-4e, 


K 6 = 


-2e. 



The null eigenvector of the reduced transition matrix M' gives the stationary probability 
distribution of equivalence classes P' st . The numbers of equivalent configurations for the 
equivalence classes ([00], [01], [02], [11], [12], [22]) are summarized as a weight vector: 

W = (l, 2, 2, 1, 2, 1) T , (A.5) 

and the stationary probability distribution of the equivalence classes P' st satisfies the nor- 
malization condition: 

W T ■ P' st = I. (A.6) 

The contributions of the classes ([00], [01], [02], [11], [12], [22]) T to the total current are de- 
scribed by 

D' = diag{0, an, a 2 , 0, 0, 0}. (A.7) 

Finally, the total current is obtained by averaging the stationary probability distribution of 
equivalence classes, weighted by the number of equivalent configurations within the classes: 

J = (W T ■ D' • P' st )/rL. (A.8) 
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